#include <vector_types.h>
#define SOLAR_BUILD_GEOMETRY_UTILS
#include <algorithm/analytical_model/algorithm/cpu/GeometryUtils.h>
#include <algorithm/analytical_model/algorithm/cuda/GeometryUtils.cuh>
#include <algorithm/analytical_model/algorithm/cuda/VectorMath.cuh>
#include <algorithm/analytical_model/algorithm/cuda/inter_op/Data.h>
#include <utils/Utils.h>

#include <cmath>

using namespace solar::cuda;
using namespace solar::cuda::inter_op;

SOLAR_CUDA_INLINE __host__ __device__ auto
solar::cuda::localToWorldRotate(const float3& d_local, const float3& aligned_normal) -> float3
{
    if (std::abs(aligned_normal.x) < float(inter_op::kEpsilon) &&
        std::abs(aligned_normal.z) < float(inter_op::kEpsilon))
    {
        return d_local;
    }
#if defined(__CUDA_ARCH__)
    float3 u, n, v, d_local1;
    n = (aligned_normal);
    u = cross(::make_float3(0.0F, 1.0F, 0.0F), n);
    u = normalize(u);
    v = cross(u, n);
    v = normalize(v);
    d_local1 = (d_local);

    float3 d_world = ::make_float3(d_local1.x * u.x + d_local1.y * n.x + d_local1.z * v.x,
                                   d_local1.x * u.y + d_local1.y * n.y + d_local1.z * v.y,
                                   d_local1.x * u.z + d_local1.y * n.z + d_local1.z * v.z);
    return d_world;
#else
    auto result = solar::cpu::localToWorldRotate(structCopy<Float3>(d_local),
                                                 structCopy<Float3>(aligned_normal));
    return structCopy<float3>(result);
#endif
}

SOLAR_CUDA_INLINE __host__ __device__ auto
solar::cuda::rotateY(const float3& origin, const float3& old_dir, const float3& new_dir) -> float3
{
#if defined(__CUDA_ARCH__)
    int dir =
        (cross(old_dir, new_dir).y > 0) ? 1 : -1; // when parallel to (0,1,0), sin = 0, no effect
    float cos = dot(old_dir, new_dir);
    float sin = dir * sqrtf(1 - cos * cos);
    return ::make_float3(cos * origin.x + sin * origin.z, origin.y,
                         -sin * origin.x + cos * origin.z);
#else
    auto result = solar::cpu::rotateY(structCopy<Float3>(origin), structCopy<Float3>(old_dir),
                                      structCopy<Float3>(new_dir));
    return structCopy<float3>(result);
#endif
}

SOLAR_CUDA_INLINE __host__ __device__ auto solar::cuda::refractRay(const float3& dir1,
                                                                   const float3& normal,
                                                                   float refract, float theta_1,
                                                                   float theta_2) -> float3
{
// r2 = (r1 / e + n1 * (cosθ1 / e – cosθ3))  *   H / cosθ3
#if defined(__CUDA_ARCH__)
    float x = (dir1.x / refract + normal.x * (std::cos(theta_1) / refract - std::cos(theta_2)));
    float y = (dir1.y / refract + normal.y * (std::cos(theta_1) / refract - std::cos(theta_2)));
    float z = (dir1.z / refract + normal.z * (std::cos(theta_1) / refract - std::cos(theta_2)));
    return ::make_float3(x, y, z);
#else
    auto result = solar::cpu::refractRay(structCopy<Float3>(dir1), structCopy<Float3>(normal),
                                         refract, theta_1, theta_2);
    return structCopy<float3>(result);
#endif
}

SOLAR_CUDA_INLINE __host__ __device__ auto
solar::cuda::rayParallelogramIntersect(const float3& orig, const float3& dir, const float3& A,
                                       const float3& B,
                                       const float3& C) -> std::tuple<float, float, float, bool>
{
#if defined(__CUDA_ARCH__)
    // #if defined(SOLAR_CUDA_PREFER_BF16)
    // bfloat16 t = NAN;
    // bfloat16 u = NAN;
    // bfloat16 v = NAN;
    // auto A_t = make_bfloat163(A);
    // auto B_t = make_bfloat163(B);
    // auto C_t = make_bfloat163(C);
    // auto dir_t = make_bfloat163(dir);
    // auto orig_t = make_bfloat163(orig);
    // auto AB = B_t - A_t;
    // auto AC = C_t - A_t;
    // // vector p
    // auto p = cross(dir_t, AC);
    // // determinant
    // auto det = dot(AB, p);

    // // if determinant is close to 0, ray lies in plane of parallelogram
    // if (__habs(det) < bfloat16(kEpsilon))
    // {
    //     return {t, u, v, false};
    // }

    // auto T = orig_t - A_t;

    // // Calculate u and make sure u <= 1
    // u = dot(T, p) / det;
    // if (u < bfloat16(0) || u > bfloat16(1))
    // {
    //     return {t, u, v, false};
    // }

    // // vector q
    // auto q = cross(T, AB);

    // // Calculate v and make sure v <= 1
    // v = dot(dir_t, q) / det;
    // if (v < bfloat16(0) || v > bfloat16(1))
    // {
    //     return {t, u, v, false};
    // }

    // t = dot(AC, q) / det;
    // return {t, u, v, t >= bfloat16(kEpsilon)};
    // #else
    float t = NAN;
    float u = NAN;
    float v = NAN;
    auto A_t = (A);
    auto B_t = (B);
    auto C_t = (C);
    auto dir_t = (dir);
    auto orig_t = (orig);
    auto AB = B_t - A_t;
    auto AC = C_t - A_t;
    // vector p
    auto p = cross(dir_t, AC);
    // determinant
    auto det = dot(AB, p);

    // if determinant is close to 0, ray lies in plane of parallelogram
    // if (abs(det) < float(kEpsilon))
    // {
    //     return {t, u, v, false};
    // }

    auto T = orig_t - A_t;

    // Calculate u and make sure u <= 1
    u = fdividef(dot(T, p), det);
    // if (u < float(0) || u > float(1))
    // {
    //     return {t, u, v, false};
    // }

    // vector q
    auto q = cross(T, AB);

    // Calculate v and make sure v <= 1
    v = fdividef(dot(dir_t, q), det);
    // if (v < float(0) || v > float(1))
    // {
    //     return {t, u, v, false};
    // }

    t = fdividef(dot(AC, q), det);
    // return {t, u, v, t >= float(kEpsilon)};
    return {t, u, v,
            std::abs(det) > float(kEpsilon) && (u >= float(0) && u <= float(1)) &&
                (v >= float(0) && v <= float(1)) && t >= float(kEpsilon)};
    // #endif
#else
    return solar::cpu::rayParallelogramIntersect(structCopy<Float3>(orig), structCopy<Float3>(dir),
                                                 structCopy<Float3>(A), structCopy<Float3>(B),
                                                 structCopy<Float3>(C));
#endif
}

SOLAR_CUDA_INLINE __host__ __device__ auto
solar::cuda::intersect(const float3& normal, const float3& plane_origin, const float3& vertex,
                       const float3& direction) -> std::tuple<float3, double>
{
#if defined(__CUDA_ARCH__)
    const auto normal_t = (normal);
    const auto direction_t1 = (direction);
    const auto plane_origin_t = (plane_origin);
    const auto vertex_t = (vertex);
    const auto div = dot(direction_t1, normal_t);
    const auto direction_t = fdividef(dot(plane_origin_t - vertex_t, normal_t), div);
    return {vertex_t + direction_t * direction_t1, direction_t};
#else
    auto result =
        solar::cpu::intersect(structCopy<Float3>(normal), structCopy<Float3>(plane_origin),
                              structCopy<Float3>(vertex), structCopy<Float3>(direction));
    auto& result_float3 = std::get<0>(result);
    return {structCopy<float3>(result_float3), std::get<1>(result)};
#endif
}
